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Abstract 


The runaway process is known to accelerate electrons in many lab- 
oratory plasmas and has been suggested as an acceleration mechanism 
in some astrophysical plasmas, including solar flares. Current calcula- 
tions of the electron velocity distributions resulting from the runaway 
process are greatly restricted because they impose spatial homogeneity 
on the distribution. We have computed runaway distributions which 
include consistent development of spatial gradients in the energetic 
tail. Our solution for the electron velocity distribution is presented as 
a function of distance along a finite length acceleration region, and is 
compared with the equivalent distribution for the infinitely long homo- 
geneous system (ie. no spatial gradients), as considered in the existing 
literature. All these results are for the weak field regime. We also 
discuss the severe restrictiveness of this weak field assumption. 



1 INTRODUCTION 


When a DC electric field is applied to a plasma runaway acceleration occurs 
of electrons above a critical velocity. This process is a consequence of the 
decrease of the dynamical friction [3] experienced by an electron when its 
speed increases above ion thermal speeds. 

This process is known to be important in laboratory plasmas such as 
tokomaks and other toroidal discharge devices [14], as well as in astrophys- 
ical plasmas eg solar figures. To date however, calculations of the velocity 
distribution and production rate of runaway electrons have been limited to 
the unphysical ideal of an homogeneous plasma of infinite spatial extent. 
This idealisation, though not strictly correct in tokomak plasmas, is rescued 
to an extent by the tokomaks toroidal periodicity. However in astrophysical 
contexts it is obviously inappropriate. 

In reality runaway acceleration occurs in finite sized acceleration regions 
which may also be inhomogeneous. These spatial gradients can significantly 
affect the development of the runaway process. In this paper we investigate 
the role of spatial gradients on the runaway process. In particular we con- 
sider how spatial gradients can develop in the tail of the distribution due to 
the finite size of the acceleration region. 

Our calculations have been made subject to the conventional assumption 
of a weak electric field. We reconsider the implications of this weak field 
assumption which is in our opinion significantly more restrictive than has 
been recognised in much of the existing literature. 

In section 2 we briefly review the existing calculations of runaway rate 
and velocity distribution. The theory and numerical technique behind our 
calculation of the electron velocity distribution is discussed in section 3. In 
section 4 we present our results and compare them with existing calculations 
which ignore the spatial variation. 


2 THEORY 

Previous attempts to understand the runaway process can be divided into 
two broad categories. The Langevin approach considers the trajectories 
of test particles moving through an unchanging background plasma. The 
second approach is to establish and solve an appropriate equation describing 
the electron velocity distribution. This second approach is clearly superior in 
principle. However it has proven so difficult in practice that it has only been 
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accomplished under the restrictive assumption that the plasma is infinite and 
homogeneous in space. 

The essence of the electron runaway process is most easily seen using the 
Langevin formulation. Consider an electron test particle moving through a 
background plasma in the presence of an externally applied DC electric field. 
The test electron is accelerated by the field. However it also experiences an 
effective frictional deceleration in the direction of its motion relative to the 
plasma. This is due to its many Coulomb interactions with the surrounding 
plasma particles. This frictional drag is greatest for an electron with velocity 
close to the thermal velocity of the local ion population [14]. Above the 
ion thermal velocity the frictional drag is a decreasing function of electron 
velocity. Thus there is always some velocity above which the electric field 
acceleration is greater than the frictional deceleration. Once above this 
critical velocity the effective acceleration of the electron grows as its velocity 
grows. This is the essence of runaway acceleration. 

The Langevin formulation is physically intuitive but incomplete. It ne- 
glects the diffusion in speed and pitch angle experienced by the test particle. 
It also lacks consistency since the background plasma can be considered as 
a collection of test particles and so should also evolve under the influence of 
the applied field. 

The electric field required to cause runaway behaviour for a thermal 
electron ( v ^ = 2kT e /m e ) is known as the Dreicer field, Ed* The importance 
of the runaway process in the overall force balance of the plasma can then 
be gauged by comparing the electric field to the Dreicer field. Calculations 
of electron distributions resulting from runaway have always been restricted 
to the weak field regime E/Ed <0.1. We also assume this restriction. This 
implies that only an exponentially small fraction of electrons will run away, 
admitting the possibility that the bulk of the distribution can remain close 
to maxwellian and near steady state. The ratio of electric field to Dreicer 
field is often rewritten as (vth/vcr) 2 , where the critical speed V&- is the speed 
at which acceleration by the electric field and frictional drag balance exactly. 

Many different estimates of the rate of production of runaway electrons 
in this weak field regime exist in the literature. The earliest were presented 
by Giovanelli [11] and Harrison [12]. However Dreicer [5, 6] is generally cred- 
ited with the first detailed analysis. He considered the abrupt application of 
a weak uniform electric field to an infinite homogeneous and ionized plasma. 
His runaway rate was determined from the rate of decay of the total popu- 
lation of collisional electrons below the critical speed. He did not determine 
the form of the distribution in the runaway regime. His solution was also 
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compromised by his assumption that the distribution function dropped to 
zero on the sphere of radius v c in velocity space. 

Both Gurevich [9] and Lebedev [17] attempted analytic solutions to the 
same problem by assuming a quasi-stationary distribution (ie. d/dt ~ 0) 
and applying different functional approximations for the electron velocity 
distribution function / over different regions of velocity space. Kruskal and 
Bernstein [15] and Gurevich and Zhivlyuk [10] developed this approach most 
fully. By identifying five seperate regions of velocity space Kruskal and Bern- 
stein determined a runaway rate to within an unknown constant multiplying 
factor. Since then a number of numerical solutions of this problem have been 
published [16, 22, 7], By comparing their numerical determination of the 
runaway rate to that of Kruskal and Bernstein, Kulsrud et al determined 
this unknown numerical multiplying factor. The resulting runaway rate, 




is the most frequently quoted runaway rate in the literature. Here u is the 
electron collision frequency and n e the electron number density. 

In the interest of clarity we point out that a critical electric field is 
often defined which is the field required to balance the drag experienced by 
an electron with speed kT e /m e . The critical and Dreicer fields are related 
by Ear = 2Ed = 47re 3 n e lnA /kT e . 


3 APPROXIMATING THE ELECTRON DISTRI- 
BUTION FUNCTION 


Our method of approximating the electron distribution function is described 
in greater detail in Ljepojevic and Burgess [18] and Ljepojevic and MacNe- 
ice [19]. The electron distribution function in a plasma can be found by 
solving the Boltzmann transport equation 


dt dx m e dv 



(1) 


where / is the electron distribution function, v is the electron velocity, x 
is the position vector, F includes all forces acting on the electrons, and the 
right hand-side of the equation describes the rate of change of / as a result 
of collisional redistribution of electrons in velocity space. It is generally true 
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that the terms on the left hand-side of equation (1) represent effects driv- 
ing the plasma away from Maxwellian, while the right hand-side describes 
thermalising collisional redistribution. 

Rosenbluth, MacDonald and Judd [20] have shown that when only 2- 
body Coulomb collisions need be considered, and v 2 /<? < 1, 

( d -L\ - _r_g_ (,?«\ , » * (, & c 

\et/ e dv» \ J dv* J 2 dv*dv* \’ dVdv" 

with 

H (v) = £ r /,(„') |r-»' r 1 . 

^ 774 ^ J 

6 J 

and 

_ 47re 4 lnA 

InA is the Coulomb logarithm, b denotes the particle species in the plasma 
(including other electrons) with which the electrons collide and c is the speed 
of light. 

Equation (1) with equation (2) is non-linear. In the low and high velocity 
regimes it admits linearisation, the Spitzer-Harm/Braginskii(SH) approxi- 
mation at low velocities, and the high velocity Landau(HVL) approximation 
at high velocities. Our solution is constructed by combining these approxi- 
mations as described below. 

For reasons of simplicity we assume that the plasma is fully ionized hy- 
drogen in steady-state. The system boundaries are two parallel planes sepa- 
rated by a distance L. Electrons entering this system across these boundaries 
are assumed to have a Maxwellian velocity distribution. A uniform electic 
field is applied between these two parallel planes in the direction perpen- 
dicular to the planes. This idealised plasma therefore will exhibit spatial 
variation only in the direction of the electric field. We shall associate a 
cartesian coordinate s with this direction. We also assume that this applied 
electric field is sufficiently weak that the lower velocity bulk of the distribu- 
tion is adequately described by the SH approximation. Finally we assume 
that there is no spatial variation in the bulk of the distribution. 
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3.1 The Spitzer-Harm/Braginskii approximation 

This approach is based on Enskog’s successive approximation (cf. Chapman 
and Cowling [4]). It proceeds under the assumptions that the distribution 
function exhibits slow temporal variation and weak spatial gradients, and is 
subject to weak electromagnetic fields. As a result the distribution function 
cam be written as 


/ — /o + /l + /2 H 7 

where successive terms represent increasingly smaller corrections. The ze- 
roth order term /o is taken to be a Maxwellian. This approximation is similar 
to a perturbation expansion in a small parameter although no small param- 
eter has been explicitly designated. In fact we shall see that the solution 
behaves as though the Knudsen parameter were that small parameter. 

Spitzer and Harm [21] considered a plasma with a weak temperature 
gradient and coaligned electric field (here assumed along the s axis). They 
adopted a first order scheme which requires 


/i = foD(v)n, 

where p = Cos# = v ■ s/v = cos# is the cosine of the electron pitch angle. D 
is a function of the electron speed determined by substituting for / in the 
Boltzmann equation, linearising the collisional integral in /i, and solving 
the resulting integro-differential equation. The solution has the form 



-Ac 


Z P E , eE 1 dp e . .Z-Dr 1 dTe] 1 

. A ' kT e + p e ds B T e day}' 


where ZDeIA and ZDt/B are tabulated in SH, p e is the electron pressure, 
T e is the electron temperature, E the electric field, A and B are normali- 
sation constants for the electric field and temperature gradient respectively, 
and Ao is the mean free path of a thermal electron. De and Dt increase 
with velocity and for p > 0 a speed v* (p) exists above which / becomes neg- 
ative. It is also clear that as the scale lengths (31nT e /ds) -1 , (dlnp e /ds)~ l 
and kT e /eE decrease, v*{p) decreases. Thus the Spitzer-Harm approach is 
only valid in the low velocity regime, and the upper boundary of this regime 
depends on the strength of the electric field and temperature and pressure 
gradients (see Gray and Kilkenny [8]). As mentioned earlier, in this paper 
we assume dT e /ds = dp e /ds = 0. 
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3.2 High velocity — the HVL approximation 

If we neglect the interaction between the high energy tail electrons in com- 
parison with their coupling to the bulk of the distribution which is close to 
Maxwellian, then we can derive an approximate linearised form of equation 
(1), applicable to the high velocity electrons. Together with the assumptions 
that the plasma is inhomogeneous in the direction of the s axis only, and in 
steady state, equation (1) then becomes [18] 


df eE df 

Vfi— - 

os m e ov 
where 


Id 2 (kT e df \1 


d_ 

dn 


i/(u) = 


47re 4 n e InA 
m\v z ’ 




( 3 ) 


with T e being the temperature characterising the bulk of the electron dis- 
tribution at height s and n e the electron number density. 

Equation (3) can then be applied to determine the velocity distribution of 
electrons at higher energies, above some appropriately chosen cutoff velocity 
v c . 

Our solution is constructed by combining the Spitzer-Harm solution be- 
low £ c = v c /vth with the HVL solution above £ c . By experimenting with 
different values £ c = 2 has been found to be the optimal value. At £ c we 
require that the HVL solution match the SH solution for nil pitch angles. 

By choosing to extend the application of the Spitzer-Harm solution out 
to £ = 2, we have placed an upper limit on the strength of the electric field 
which we can consider. We will choose to set E/E D = 0.05. At this field 
strength the Spitzer-Harm distribution becomes negative for 0 = n and f in 
excess of about 2.2. At higher field strengths the Spitzer-Harm distribution 
becomes negative at even lower speeds. 

Equation (3) is a linearised form of the non-relativistic Landau equation. 
Obviously when the electron velocity becomes a significant fraction of the 
speed of light relativistic effects become important and equation (3) becomes 
invalid. An appropriate relativistic generalisation of the Landau equation 
has been derived [1, 13]. However it is considerably more complex in form, 
and for reasons of simplicity we have chosen to retain the non-relativistic 
form for this present work. This clearly limits the maximum temperatures 
of the bulk distributions which we can consider. 
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Figure 1: ln(///j\r) variation with speed ( and pitch angle 9 at s = 4.375 x 
10 9 cm, for 0 < 9 < n and 0 < £ < 25. Unlabelled contours are evenly spaced 
between labelled contours. 

4 RESULTS 

Equation (3) was solved for a plasma with bulk temperature 10 6 K and elec- 
tron number density 10 9 cm~ 3 . A uniform electric field E = — 0.05i?£)S 
was applied between the plane parallel system boundaries separated by a 
distance L = 5 x 10 9 cm. The solution was obtained for electron speeds 
up to £ = 25 on a finite difference mesh with 80 grid cells along the 
s axis. The mean free path of a thermal electron in this plasma is ap- 
proximately 8 x 10 6 cm, or approximately one tenth of the spatial grid cell 
size. Since the electron mean free path scales as f 4 , there are roughly two 
grid cells per mean free path at £ = 2. The normalised critical velocity is 
&r = W”t/» = (£W^ 1/2 ~4.5. 
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Figure 2: (a) Variation of ln(f/f N ) with speed £ for pitch angle 9 = 0, ie. 
in the direction of electron acceleration, at s = 4.375 x 10 9 cm. The dashed 
curve is the equivalent maxwellian curve, (b) The same data plotted in (a) 
recast in the form of the exponent function g. 

4.1 The Distribution Function 

A convenient representation of the solution can be achieved by defining a 
function p(s,£,/x) such that 

/ = f N ex p (-pC 2 ) 

where the normalisation constant f N = n{2irkT e /m e )~ 3 / 2 . A value of g = 
1 gives the maxwellian value for f,g> 1 underpopulation relative to a 
maxwellian and g < 1 overpopulation relative to a maxw ellian 

Figure 1 illustrates the structure in velocity space of our solution at a 
selected location, s = 4.375 x 10 9 cm. The tail of the distribution is over- 
populated relative to a Maxwellian distribution for 6 <n/2. The variation 
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of both ln(///Ar) and g with £ for 9 = 0, ie the direction anti-parallel to 
the applied field, is shown in figures 2 a and b. At low velocity and very 
high velocity the distribution is close to maxwellian with temperature equal 
to the assumed bulk temperature. The distribution begins to depart sig- 
nificantly from maxwellian near £ ~ 3.5 as the increasingly collisionless 
character of the electrons permits them to gain more energy from the ap- 
plied field between collisions. The ln(///w) curve “flattens” out somewhat 
before eventually turning down again to rejoin the original maxwellian curve 
at very high velocities. The extent of this “flattened” section of the curve is 
limited by the finite spatial extent of the acceleration region. As the length 
of the acceleration region is extended this flattened section is also extended 
(figure 3) and would ultimately extend to infinite velocities in the idealised 
case of an infinitely long acceleration region. 

The typical variation of In/ with 6 in the tail is illustrated in figure 4. 
Note that the distribution is enhanced above maxwellian in the direction 
perpendicular to the field. This is a consequence of the diffusive spread in 
pitch angle of the runaway tail. 

In figure 5 we have plotted the variation of In/ with perpendicular speed 
fx at s = 4.375 x 10 9 cm for two different values of = £ - $ of 6 and 10. 
For < 5 where most of the energy resides, this is well fit by a gaussian 
with a temperature of 1.6 x 10 6 K, as indicated by the superimposed dashed 
line in figure 5b. Broadenning of this perpendicular velocity distribution 
persists at larger £x but with steadily diminishing degree, until for £x > 20 
the perpendicular distribution becomes Maxwellian at about the original 
background temperature of 10 6 K. To determine the energy content of the 
perpendicular velocity distribution we define a perpendicular temperature 

t±((\\,s) = J~ dt ± e±m,u,s)/ J~ dt ± um\\,u,s) 

following Fuchs et al [7]. This is plotted in figure 6 as a function of at a 
number of different spatial locations throughout the acceleration region. For 
example at s = 4.375 x 10 9 cm and = 10 the perpendicular temperature 
is 1.6 x 10 6 K as we would expect given the fit shown in figure 5b. 

Obviously, in the absence of any explicit source of electrical current there 
should be no spatial variation of the total current. In our solution however 
the spatial variation of / above £ = 2 represents a small variation in current. 
In reality any variation of this kind in the current would produce a local 
accumulation of charge resulting in a local electric field acting to disperse 
this accumulation. Therefore, in applying a fixed potential between the 
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Figure 3: Vaxiation of In (///at) with speed £ for pitch angle 0 = 0 beginning 
at s = Ocm (leftmost solid curve), and continuing at regular intervals of 
6.25 x 10 8 cm up to 4.375 x 10 9 cm (rightmost solid curve). The dotted line 
is for the infinitely long homogeneous case. 


boundaries of the system we have not guaranteed that a uniform field will 
exist across the system. However in practice these departures from constant 
current which we observe are so slight that the local fields which should 
develop to restore current conservation can be neglected in comparison with 
the much larger applied field. The variation of the current density with 
position is illustrated in table 1. 

The current flowing at s 0 = 0 is -3.364 x 10 7 statamps cm' 2 , while 
the current at s = 5 x 10 9 cm is -3.436 x 10 7 statamps cm' 2 . Most of the 
current is carred by the electrons below £ = 2 which contribute —2.514 x 
10 7 statamps cm 2 to the total current density. This enables us to use the 
Spitzer-Harm electrical conductivity to calculate the approximate field cor- 
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Figure 4: Variation of In {f / /n) with pitch angle 6 for velocity £ = 7s at 
s = 6.25 x 10 8 cm (lower curve), s = 2.5 x 10 9 cm and s = 4.375 x 10 9 cm 
(upper curve). 


rection necessary to ensure constant current. In fact the ratio AE/ E{sq) 
should be approximately equal to Aj/j(so) = 0.02. This correction would 
have a negligible effect on the tail of the distribution. 

What difference does inclusion of spatial gradients make to the velocity 
distribution? To answer this question we need to compare our results with 
the solution obtained when d/ds = 0, ie. the system is assumed infinite and 
homogeneous. As we have mentioned, such calculations have already been 
published [22, 16, 7]. Unfortunately these authors have always considered 
E/Ed > 0.08, and have applied the same linearised equation throughout 
all of velocity space with no special consideration of the thermal velocity 
regime. We have been limited to E/Ed < 0.05 precisely because of our 
application of the Spitzer-Harm distribution at thermal velocities. Therefore 
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Figure 5: Distribution of perpendicular velocities of electrons at s = 4.375 x 
10 9 cm (solid lines) with (a) fy = £ • s = 6 and (b) = 10. The dotted 

curves are the equivalent result for the infinite homogeneous case. 


to isolate the differences introduced by consideration of spatial gradients it 
was necessary for us to generate a solution for the infinite homogeneous 
distribution with E/Ed — 0.05 and consistent with the application of the 
Spit zer- Harm distribution below £ = 2. We used a relaxation approach 
similar to that of Wiley at al. [22]. The appropiate equation, 


?l-?E.dL-± d 2 , x/*r«a/ \i a r /t 2 ,an 

dt m e dv v 2 dv [ U \m e dv Vj dp [ ^ 5/i ’ 

was solved using a fimte difference approximation until an almost steady 
state was achieved. This steady state distribution is illustrated in figures 3, 
5 and 6. The infinite homogeneous solution appears to be the asymptotic 
limit to which the position dependent distribution tends as the length of the 
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Figure 6: Variation of the perpendicular temperature T± with £|| and spatial 
location. Each curve represents a different spatial location beginning with 
the lowest solid curve at s = 6.25 x 10 8 cm, and continuing at regular intervals 
of 6.25 x 10 8 cm out to the highest solid curve at s = 4.375 x 10 9 cm. The 
dotted curve is for the infinitely long homogeneous case. 

acceleration region approaches infinity. The tail of this limi t distribution, 
for 0 = 0, has an / oc l/£ dependence, as expected from conservation of flux 
arguments [2]. 

4.2 The Runaway Rate 

How do these results compare with the runaway rates quoted in the litera- 
ture? To compute runaway rates we have followed the example of Wiley et 
al [22]. Let us consider first the infinite homogeneous system. The rate of 
production of runaways can be computed from the rate of loss of electrons 
from the collisional bulk of the distribution. We choose to idealise the colli- 
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Table 1: Electric Current Density as a function of position for a uniform 
electric field E = —0.05 Eds. 


s (10 9 cm) 

j (10 7 statamps.cm 2 ) 

0.0000 

3.364 

0.1250 

3.421 

0.2500 

3.422 

0.3750 

3.417 

0.5000 

3.414 

0.6250 

3.413 

0.7500 

3.413 

4.6250 

3.413 

4.7500 

3.413 

4.8750 

3.416 

5.0000 

3.436 


sional regime as a sphere in velocity of radius £r. Integrating the equation 
(5) over the volume of the collisional regime gives 


dnc dun 87r 2 ne 4 lnA 


dt 


at 


mi 


!'-M 


E 

Ed 


Ckfi^R.)^ + 


1 df 


2 ZiidZ 


+ 


R 


/(£*)) 


( 6 ) 


where nc is the number of electrons in the bulk and hr is the number of 
runaway electrons. For comparative purposes we define a normalised time 


47re 4 lnA 1 

T==< 2 7=- 

m% 7T0T 


The rimaway rate T = n~ l dnR/dr can be easily calculated once we have 
chosen an appropriate value for The obvious choice would seem to be 
the critical velocity The variation of T with ^ is illustrated in 

figure 7. For the infinite homogeneous system, the rate decreases rapidly 
with increasing until just above the critical velocity, near 4.5 in this case. 
It maintains a fairly constant value thereafter. This constant rate, which 
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Figure 7: Variation of the rate T as a function of for the case of 
the infinitely long homogeneous plasma(dotted line), and for position s = 
4.375 x 10 9 cm along the finite sized acceleration region(solid line). The hori- 
zontal broken line indicates the runaway rate computed from the expression 
of Kruskal and Bernstein given in section 2. 


is approximately 70% larger than predicted by the normalised Kruskal and 
Bernstein [15] formula is the obvious choice for the runaway rate. 

When spatial gradients are retained in the calculation of the distribution 
function the identification of a runaway rate is much less clear. The integral 
of equation (3) over the collisional sphere in velocity space produces an 
equation identical to equation (6) except for the addition of a term due to 
the spatial gradient. However the spatial gradient term is responsible for 
convective transport of electrons, and does not result in any change to their 
speed. It does not contribute directly to the rate of production of runaways. 
Thus equation (6) can be used to compute the runaway rate in the presence 
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Position ( x 10 9 cm) 


Figure 8: Variation of the runaway rate with position in the finite length 
acceleration region. The rate is computed using = 5. 


of spatial gradients also. In figure 7 we show the variation of T with £r at 
s — 4.375 x 10 9 cm. As £r increases the rate drops steadily. The energetic 
tail has not developed sufficiently to produce a rate independent of £r as was 
the case in the infinite homogeneous system. The runaway rate is clearly 
sensitively dependent on our choice of £r. We will choose £r = 5 for two 
reasons. The first is that it lies just above the critical runaway velocity 
£cr = 4.5 which is where the rate became insensitive to the value of £r in 
the infinite homogeneous case. The second is that this is the same choice 
adopted by Wiley et al [22]. In figure 8 we show how this runaway rate 
varies as a function of position along the acceleration region. 
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5 DISCUSSION 

This present calculation is limited in the maximum values of E / Eq and T e 
which can be considered. Our requirement to match solutions above and 
below £ = 2 limits the size of the electric field . For E/Ed > 05 the 
Spitzer-Harm solution is no longer sufficiently accurate near £ = 2. T his 
limitation can only be overcome by more accurately solving the non-linear 
equation which defines the distribution function in the region of the thermal 
velocity. No other numerical calculations of the runaway rate have addressed 
this difficulty satisfactorily in spite of the fact that they typically apply 
electric fields in the range 0.08 < E/E D < 0.2. The limitation on the bulk 
temperatures which we can consider is a consequence of our non-relativistic 
treatment. By raising the temperature we push the runaway regime to a 
higher range in (u/c) 2 . 

Our results show that except during the initial stages of acceleration 
the runaway rate predicted by the normalised Kruskal and Bernstein theory 
is accurate to within a factor of two. However in any application which 
requires the actual distribution of electrons rather than the total number 
in the tail, the infinite homogeneous result is capable of being in error by 
orders of magnitude. 
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